rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#Make sure that your working folder has two sub-folders:
#-data
#-paper

#With paper, you should also have two sub-folders:
#-tables
#-graphs

#Please install any of the libraries listed below that are not installed on your machine.

#This is to obtain:
#-table_1
#-table_A1
#-figure_a12a
#-figure_a12b
#-figure_a12c

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("fastDummies")
library("dplyr")
library("berryFunctions")
library("mediation")
library("sp")
library("pgirmess")
library("fwildclusterboot")
library("geojsonsf")
library("tidyverse")

#Functions

dlshape=function(shploc, shpfile) {
  temp=tempfile()
  download.file(shploc, temp)
  unzip(temp)
  fp <- file.path(temp)
  obj_shp<-rgdal::readOGR(".",shpfile)
  return(obj_shp)}

mean_sd_fun<-function(lm_model) {
  list_empty <-c()
  number_mean<-mean(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  number_sd<-sd(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  list_empty[1]<-number_mean
  list_empty[2]<-number_sd
  names(list_empty)<-c("mean", "sd")
  return(list_empty)
}


week_obs_count<-function(lm_model) {
  list_empty <-c()
  week_no<-length(grep(x = colnames(data.frame(lm_model$model)), pattern = "^week"))
  #Adding 2 because there should be 53 weeks and because I am taking 2 weeks out as in STATA 2 weeks drop
  week_no<-week_no+2
  total_obs<-nobs(lm_model)
  #counties<-ncol(lm_model$model%>% dplyr:: select(starts_with("AGS_")))
  counties<-total_obs/week_no
  list_empty[1]<-week_no
  list_empty[2]<-total_obs
  list_empty[3]<-counties
  names(list_empty)<-c("week_no", "total_obs", "counties")
  return(list_empty)
}


#Read CSV file
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")

#############################
#Function for calculating SE#
#############################

#Step1. Turning some variables to numeric
data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)

data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
summary(data_cov_mob$bridging_tot_pop)
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)


#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
summary(data_cov_mob$bonding_tot_pop)
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)

#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000
summary(data_cov_mob$vereine_tot_pop)
mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)

#Right Wing
data_cov_mob$pct_afd<-as.numeric(data_cov_mob$pct_afd)
summary(data_cov_mob$pct_afd)
mean_pct_afd<-summary(data_cov_mob$pct_afd)[3]
data_cov_mob$pct_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_afd, 1, 0)

#Right Wing Change
data_cov_mob$pct_change_afd<-as.numeric(data_cov_mob$pct_change_afd)
summary(data_cov_mob$pct_change_afd)
mean_pct_change_afd<-summary(data_cov_mob$pct_change_afd)[4]
data_cov_mob$pct_change_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_change_afd, 1, 0)


table(data_cov_mob$vereine_tot_pop_dummy)
table(data_cov_mob$bridging_tot_pop_dummy)
table(data_cov_mob$pct_afd_dummy)
table(data_cov_mob$pct_change_afd_dummy)


data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop", 
                                                        "bonding_tot_pop",
                                                        "pct_afd_dummy",
                                                        "pct_change_afd_dummy"))
summary(data_cov_mob$gender_ratio)
data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)

data$post_treat_vereine_tot_pop_dummy<-data$post_treat*data$vereine_tot_pop_dummy
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy
data$post_treat_bridging_tot_pop<-data$post_treat*data$bridging_tot_pop

data$post_treat_bonding_tot_pop_dummy<-data$post_treat*data$bonding_tot_pop_dummy
data$post_treat_pct_afd_dummy<-data$post_treat*data$pct_afd_dummy
data$post_treat_pct_change_afd_dummy<-data$post_treat*data$pct_change_afd_dummy
data$post_treat_pct_change_afd<-data$post_treat*data$pct_change_afd

unique(data$post_treat_vereine_tot_pop_dummy)
unique(data$post_treat_bridging_tot_pop_dummy)
unique(data$post_treat_bonding_tot_pop_dummy)
unique(data$post_treat_pct_afd_dummy)
unique(data$post_treat_pct_change_afd_dummy)


#Save to separate file
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(nuts2$AGS)

# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')
summary(datas$gender_ratio)


new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))


list_dummies<-names(new)
#I am removing the counties which get dropped in Stata
#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
list_week_dummies_now1<-list_week_dummies[list_week_dummies!= "week_01"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_10"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_53"]


list_land_dummies<-names(new3)
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]
temp_lag<-10

#Creating the interaction terms for bonding
for (i in list_land_dummies){
  for (j in list_week_dummies_now1)
    datas[[paste("inter_", i, "_", j, sep="")]]<-datas[[i]]*datas[[j]]
}

new<-datas %>% dplyr::select(starts_with("inter_"))
list_land_dummies_star<-names(new)

#VEREINE ALL
basic_ver<-c("vereine_tot_pop_dummy", "post_treat", "post_treat_vereine_tot_pop_dummy", "lag_covid_pc")

extended_ver <- c("vereine_tot_pop_dummy", "post_treat", "post_treat_vereine_tot_pop_dummy", "lag_covid_pc",
                  "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                  "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                  "east_ger", "gender_ratio", "pct_students")


basic_and_dummies_ver<-list.append(basic_ver, list_dummies, list_week_dummies2)
basic_and_dummies_int_ver<-list.append(basic_ver, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_ver<-list.append(extended_ver, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_ver<-list.append(extended_ver, list_land_dummies, list_week_dummies2, list_land_dummies_star)


spec_basic_ver<-as.formula(paste("mean_mobil~",
                                 paste(basic_and_dummies_ver, collapse= "+")))
spec_basic_int_ver<-as.formula(paste("mean_mobil~",
                                     paste(basic_and_dummies_int_ver, collapse= "+")))

spec_extended_ver<-as.formula(paste("mean_mobil~",
                                    paste(extended_and_week_dummies_ver, collapse= "+")))
spec_extended_int_ver<-as.formula(paste("mean_mobil~",
                                        paste(extended_and_week_dummies_int_ver, collapse= "+")))


#BRIDGING
basic_brid<-c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc")
extended_brid <- c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students")

basic_and_dummies_brid<-list.append(basic_brid, list_dummies, list_week_dummies2)
basic_and_dummies_int_brid<-list.append(basic_brid, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2, list_land_dummies_star)

spec_basic_brid<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_brid, collapse= "+")))
spec_basic_int_brid<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_brid, collapse= "+")))

spec_extended_brid<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_brid, collapse= "+")))
spec_extended_int_brid<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_brid, collapse= "+")))

#BONDING
basic_bond<-c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc")
extended_bond <- c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students")

basic_and_dummies_bond<-list.append(basic_bond, list_dummies, list_week_dummies2)
basic_and_dummies_int_bond<-list.append(basic_bond, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_bond<-list.append(extended_bond, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_bond<-list.append(extended_bond, list_land_dummies, list_week_dummies2, list_land_dummies_star)


spec_basic_bond<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_bond, collapse= "+")))
spec_basic_int_bond<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_bond, collapse= "+")))

spec_extended_bond<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_bond, collapse= "+")))
spec_extended_int_bond<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_bond, collapse= "+")))



temp_lag<-10

#*Model1 - Stata
#xtreg mean_mobil vereine_tot_pop_dummy /*
#  */post_treat /*
#  */post_treat_vereine_tot_pop_dummy /*
#  */lag_covid_pc /*
#  */ags_* /*
#  */week_* /*
#  */i.sn_l##i.week
#Note the variables in Stata which get dropped

datas_city<-subset(datas, BEZ=="Kreisfreie Stadt")
datas_nocity<-subset(datas, BEZ=="Kreis")
mean(datas$pop_den, na.rm=T)
datas_above_meandens<-subset(datas, pop_den>mean(datas$pop_den, na.rm=T))
datas_below_meandens<-subset(datas, pop_den<mean(datas$pop_den, na.rm=T))


mo1_verx<-lm(spec_basic_int_ver, data=datas_city)
summary(mo1_verx)

mo1_ver<-lm(spec_basic_int_ver, data=datas)
summary(mo1_ver)
mo1_ver_df<-broom::tidy(mo1_ver)
mo1_ver_df$rse<-coeftest(mo1_ver, vcov = cluster.vcov(mo1_ver, datas$AGS, force_posdef=T))[, 2]
mo1_ver_df$pv_rse<-coeftest(mo1_ver, vcov = cluster.vcov(mo1_ver, datas$AGS, force_posdef=T))[, 4]
mo1_ver_df2<-mo1_ver_df %>% filter(term %in% basic_ver)

x3<-labels(terms(spec_basic_int_ver))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$gvereine
names(datas3_oneweek)


#Conley SE
#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts$AGS<-as.numeric(nuts$AGS)
nuts_merged<-left_join(nuts, datas3_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$gvereine))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
nuts_merged3$gvereine<-as.numeric(nuts_merged3$gvereine)
pgi_cor <- correlog(coords=coo, z=nuts_merged3$gvereine, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
pgi_cor_df$dist.class


#Saving the correlogram
moran_gvereine<-ggplot(data=pgi_cor_df, aes(x=dist.class/100000, y=coef, group=1)) +
  geom_line()+
  geom_point()+
  geom_hline(yintercept = 0, color ="red")+
  labs(x = "Distance Class", y = "Moran I statistic")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(moran_gvereine, file = "paper/graphs/figure_a12a.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

#Calculating the distance the minimises Moran's I stat
dist_cut_ver<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_ver, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)

res_basic_int_ver<-conleyreg(spec_basic_int_ver, data=datas3, dist_cutoff=dist_cut_ver, unit = "AGS", time = "time", 
                             dist_mat=dm_bri, lag_cutoff=temp_lag, rowwise=T)

res_basic_int_ver_df<-broom::tidy(res_basic_int_ver)
res_basic_int_ver_df2<-res_basic_int_ver_df %>% filter(term %in% basic_ver)
res_basic_int_ver_df2<-subset(res_basic_int_ver_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_ver_df2)<-c("term", "con_std_error", 'con_p_value')

model1<-left_join(mo1_ver_df2, res_basic_int_ver_df2, by = c("term" = "term"))

#Creating a folder to save table1
dir.create(file.path(working_folder, "./paper/tables/table1/"), showWarnings = FALSE)
write.csv(model1,"./paper/tables/table1/model1_ver.csv", row.names = FALSE)

#Model2 - Stata
#Very close but slightly different
#xtreg mean_mobil vereine_tot_pop_dummy /*
#  */post_treat /*
#  */post_treat_vereine_tot_pop_dummy /*
#  */lag_covid_pc /*
#  *ags_*
#  */week_* /*
#  */i.sn_l##i.week /*
# */ pct_turnout log_pop_total log_gdp_per_cap pop_den /*
#  */ pct_college pct_pop_above_60 pct_service_sec_narrow pct_manufacturing



mo8_ver<-lm(spec_extended_int_ver, data=datas)
summary(mo8_ver)
mo8_ver_df<-broom::tidy(mo8_ver)
mo8_ver_df$rse<-coeftest(mo8_ver, vcov = cluster.vcov(mo8_ver, datas$AGS, force_posdef=T))[, 2]
mo8_ver_df$pv_rse<-coeftest(mo8_ver, vcov = cluster.vcov(mo8_ver, datas$AGS, force_posdef=T))[, 4]
mo8_ver_df2<-mo8_ver_df %>% filter(term %in% extended_ver)


x3<-labels(terms(spec_extended_int_ver))
datas$time<-as.numeric(datas$week)
datas3<-subset(datas, select = list.append(x3, "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]


dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_ver, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_extended_int_ver, data=datas3, dist_cutoff=dist_cut_ver, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_extended_int_ver_df<-broom::tidy(res)
res_extended_int_ver_df2<-res_extended_int_ver_df %>% filter(term %in% extended_ver)
res_extended_int_ver_df2<-subset(res_extended_int_ver_df2, select = c("term", "std.error", 'p.value'))
names(res_extended_int_ver_df2)<-c("term", "con_std_error", 'con_p_value')

model2<-left_join(mo8_ver_df2, res_extended_int_ver_df2, by = c("term" = "term"))
write.csv(model2,"./paper/tables/table1/model2_ver.csv", row.names = FALSE)


#*Model3
#*County FE
#xtreg mean_mobil bridging_tot_pop_dummy /*
#  */post_treat /*
#  */post_treat_bridging_tot_pop_d /*
#  */lag_covid_pc /*
#  */ags_* /*
#  */week_* /*
#  */i.sn_l##i.week


mo1_brid<-lm(spec_basic_int_brid, data=datas)
summary(mo1_brid)
mo1_brid_df<-broom::tidy(mo1_brid)
mo1_brid_df$rse<-coeftest(mo1_brid, vcov = cluster.vcov(mo1_brid, datas$AGS, force_posdef=T))[, 2]
mo1_brid_df$pv_rse<-coeftest(mo1_brid, vcov = cluster.vcov(mo1_brid, datas$AGS, force_posdef=T))[, 4]
mo1_brid_df2<-mo1_brid_df %>% filter(term %in% basic_brid)

basic_and_dummies_int_brid_x<-list.append(basic_brid, list_land_dummies_star)
spec_basic_int_brid_x<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_brid_x, collapse= "+"), "| AGS + week"))

summary(mo1_brid)


x3<-labels(terms(spec_basic_int_brid))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)

#Saving the correlogram
moran_bridging_total<-ggplot(data=pgi_cor_df, aes(x=dist.class/100000, y=coef, group=1)) +
  geom_line()+
  geom_point()+
  geom_hline(yintercept = 0, color ="red")+
  labs(x = "Distance Class", y = "Moran I statistic")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(moran_bridging_total, file = "paper/graphs/figure_a12b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_brid, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_brid<-broom::tidy(res)
res_basic_int_brid_df2<-res_basic_int_brid %>% filter(term %in% basic_brid)
res_basic_int_brid_df2<-subset(res_basic_int_brid_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model3<-left_join(mo1_brid_df2, res_basic_int_brid_df2, by = c("term" = "term"))
write.csv(model3,"./paper/tables/table1/model3_brid.csv", row.names = FALSE)


# *Model4
# *NO County FE
# *Week FE
# xtreg mean_mobil bridging_tot_pop_dummy /*
#   */post_treat /*
#   */post_treat_bridging_tot_pop_d /*
#   */lag_covid_pc /*
#   *ags_*
#   */week_* /*
#   */i.sn_l##i.week /*
# */ pct_turnout log_pop_total log_gdp_per_cap pop_den /*
#   */ pct_college pct_pop_above_60 pct_service_sec_narrow pct_manufacturing


mo8_brid<-lm(spec_extended_int_brid, data=datas)
summary(mo8_brid)

mo8_brid_df<-broom::tidy(mo8_brid)
mo8_brid_df$rse<-coeftest(mo8_brid, vcov = cluster.vcov(mo8_brid, datas$AGS, force_posdef=T))[, 2]
mo8_brid_df$pv_rse<-coeftest(mo8_brid, vcov = cluster.vcov(mo8_brid, datas$AGS, force_posdef=T))[, 4]
mo8_brid_df2<-mo8_brid_df %>% filter(term %in% extended_brid)


x3<-labels(terms(spec_extended_int_brid))
datas$time<-as.numeric(datas$week)

datas3<-subset(datas, select = list.append(x3, "mean_mobil", "bridging_total", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_extended_int_brid, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_extended_int_brid_df<-broom::tidy(res)
res_extended_int_brid_df2<-res_extended_int_brid_df %>% filter(term %in% extended_brid)
res_extended_int_brid_df2<-subset(res_extended_int_brid_df, select = c("term", "std.error", 'p.value'))
names(res_extended_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model4<-left_join(mo8_brid_df2, res_extended_int_brid_df2, by = c("term" = "term"))
write.csv(model4,"./paper/tables/table1/model4_brid.csv", row.names = FALSE)

# *Model5
# *County FE
# xtreg mean_mobil bonding_tot_pop_dummy /*
#   */post_treat /*
#   */post_treat_bonding_tot_pop_d /*
#   */lag_covid_pc /*
#   */ags_* /*
#   */week_* /*
#   */i.sn_l##i.week


mo1_bond<-lm(spec_basic_int_bond, data=datas)
summary(mo1_bond)
mo1_bond_df<-broom::tidy(mo1_bond)
mo1_bond_df$rse<-coeftest(mo1_bond, vcov = cluster.vcov(mo1_bond, datas$AGS, force_posdef=T))[, 2]
mo1_bond_df$pv_rse<-coeftest(mo1_bond, vcov = cluster.vcov(mo1_bond, datas$AGS, force_posdef=T))[, 4]
mo1_bond_df2<-mo1_bond_df %>% filter(term %in% basic_bond)


x3<-labels(terms(spec_basic_int_bond))
datas$time<-as.numeric(datas$week)


datas3<-subset(datas, select = list.append(x3, "bonding_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bonding_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bonding_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)

#Saving the correlogram
moran_bonding_total<-ggplot(data=pgi_cor_df, aes(x=dist.class/100000, y=coef, group=1)) +
  geom_line()+
  geom_point()+
  geom_hline(yintercept = 0, color ="red")+
  labs(x = "Distance Class", y = "Moran I statistic")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(moran_bonding_total, file = "paper/graphs/figure_a12c.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



#Calculating the distance the minimises Moran's I stat
dist_cut_bon<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_bon, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_bond, data=datas3, dist_cutoff=dist_cut_bon, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_bond_df<-broom::tidy(res)
res_basic_int_bond_df2<-res_basic_int_bond_df %>% filter(term %in% basic_bond)
res_basic_int_bond_df2<-subset(res_basic_int_bond_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_bond_df2)<-c("term", "con_std_error", 'con_p_value')

model5<-left_join(mo1_bond_df2, res_basic_int_bond_df2, by = c("term" = "term"))
write.csv(model5,"./paper/tables/table1/model5_bon.csv", row.names = FALSE)

# *Model6
# *NO County FE
# *Week FE
# xtreg mean_mobil bonding_tot_pop_dummy /*
#   */post_treat /*
#   */post_treat_bonding_tot_pop_d /*
#   */lag_covid_pc /*
#   *ags_*
#   */week_* /*
#   */i.sn_l##i.week /*
# */ pct_turnout log_pop_total log_gdp_per_cap pop_den /*
#   */ pct_college pct_pop_above_60 pct_service_sec_narrow pct_manufacturing

mo8_bond<-lm(spec_extended_int_bond, data=datas)
summary(mo8_bond)
mo8_bond_df<-broom::tidy(mo8_bond)
mo8_bond_df$rse<-coeftest(mo8_bond, vcov = cluster.vcov(mo8_bond, datas$AGS, force_posdef=T))[, 2]
mo8_bond_df$pv_rse<-coeftest(mo8_bond, vcov = cluster.vcov(mo8_bond, datas$AGS, force_posdef=T))[, 4]
mo8_bond_df2<-mo8_bond_df %>% filter(term %in% extended_bond)

x3<-labels(terms(spec_extended_int_bond))
datas$time<-as.numeric(datas$week)


datas3<-subset(datas, select = list.append(x3, "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]

temp_lag<-10

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_bon, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_extended_int_bond, data=datas3, dist_cutoff=dist_cut_bon, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_extended_int_bond_df<-broom::tidy(res)
res_extended_int_bond_df2<-res_extended_int_bond_df %>% filter(term %in% extended_bond)
res_extended_int_bond_df2<-subset(res_extended_int_bond_df2, select = c("term", "std.error", 'p.value'))
names(res_extended_int_bond_df2)<-c("term", "con_std_error", 'con_p_value')

model6<-left_join(mo8_bond_df2, res_extended_int_bond_df2, by = c("term" = "term"))
write.csv(model6,"./paper/tables/table1/model6_bon.csv", row.names = FALSE)

#Organizing Table

mo1_ver<-lm(spec_basic_int_ver, data=datas)
mo8_ver<-lm(spec_extended_int_ver, data=datas)
mo1_brid<-lm(spec_basic_int_brid, data=datas)
mo8_brid<-lm(spec_extended_int_brid, data=datas)
mo1_bond<-lm(spec_basic_int_bond, data=datas)
mo8_bond<-lm(spec_extended_int_bond, data=datas)


model1_csv <- as.data.frame(read_csv("./paper/tables/table1/model1_ver.csv"))
model2_csv <- as.data.frame(read_csv("./paper/tables/table1/model2_ver.csv"))
model3_csv <- as.data.frame(read_csv("./paper/tables/table1/model3_brid.csv"))
model4_csv <- as.data.frame(read_csv("./paper/tables/table1/model4_brid.csv"))
model5_csv <- as.data.frame(read_csv("./paper/tables/table1/model5_bon.csv"))
model6_csv <- as.data.frame(read_csv("./paper/tables/table1/model6_bon.csv"))

model1<-model1_csv
model2<-model2_csv
model3<-model3_csv
model4<-model4_csv
model5<-model5_csv
model6<-model6_csv


model1$order<-NA
model2$order<-NA
model3$order<-NA
model4$order<-NA
model5$order<-NA
model6$order<-NA

model1$order[model1$term=="vereine_tot_pop_dummy"]<-1
model2$order[model2$term=="vereine_tot_pop_dummy"]<-1
model3$order[model3$term=="bridging_tot_pop_dummy"]<-2
model4$order[model4$term=="bridging_tot_pop_dummy"]<-2
model5$order[model5$term=="bonding_tot_pop_dummy"]<-3
model6$order[model6$term=="bonding_tot_pop_dummy"]<-3

model1$order[model1$term=="post_treat"]<-4
model2$order[model2$term=="post_treat"]<-4
model3$order[model3$term=="post_treat"]<-4
model4$order[model4$term=="post_treat"]<-4
model5$order[model5$term=="post_treat"]<-4
model6$order[model6$term=="post_treat"]<-4

model1$order[model1$term=="post_treat_vereine_tot_pop_dummy"]<-5
model2$order[model2$term=="post_treat_vereine_tot_pop_dummy"]<-5
model3$order[model3$term=="post_treat_bridging_tot_pop_dummy"]<-6
model4$order[model4$term=="post_treat_bridging_tot_pop_dummy"]<-6
model5$order[model5$term=="post_treat_bonding_tot_pop_dummy"]<-7
model6$order[model6$term=="post_treat_bonding_tot_pop_dummy"]<-7


model2$order[model2$term=="pct_turnout"]<-8
model4$order[model4$term=="pct_turnout"]<-8
model6$order[model6$term=="pct_turnout"]<-8

model2$order[model2$term=="log_pop_total"]<-9
model4$order[model4$term=="log_pop_total"]<-9
model6$order[model6$term=="log_pop_total"]<-9

model2$order[model2$term=="log_gdp_per_cap"]<-10
model4$order[model4$term=="log_gdp_per_cap"]<-10
model6$order[model6$term=="log_gdp_per_cap"]<-10

model2$order[model2$term=="pop_den"]<-11
model4$order[model4$term=="pop_den"]<-11
model6$order[model6$term=="pop_den"]<-11

model2$order[model2$term=="pct_college"]<-12
model4$order[model4$term=="pct_college"]<-12
model6$order[model6$term=="pct_college"]<-12

model2$order[model2$term=="pct_pop_above_60"]<-13
model4$order[model4$term=="pct_pop_above_60"]<-13
model6$order[model6$term=="pct_pop_above_60"]<-13

model2$order[model2$term=="pct_pop_under_35"]<-14
model4$order[model4$term=="pct_pop_under_35"]<-14
model6$order[model6$term=="pct_pop_under_35"]<-14

model2$order[model2$term=="pct_service_sec_narrow"]<-15
model4$order[model4$term=="pct_service_sec_narrow"]<-15
model6$order[model6$term=="pct_service_sec_narrow"]<-15

model2$order[model2$term=="pct_manufacturing"]<-16
model4$order[model4$term=="pct_manufacturing"]<-16
model6$order[model6$term=="pct_manufacturing"]<-16

model2$order[model2$term=="east_ger"]<-17
model4$order[model4$term=="east_ger"]<-17
model6$order[model6$term=="east_ger"]<-17

model1$order[model1$term=="lag_covid_pc"]<-18
model2$order[model2$term=="lag_covid_pc"]<-18
model3$order[model3$term=="lag_covid_pc"]<-18
model4$order[model4$term=="lag_covid_pc"]<-18
model5$order[model5$term=="lag_covid_pc"]<-18
model6$order[model6$term=="lag_covid_pc"]<-18


model2$order[model2$term=="gender_ratio"]<-19
model4$order[model4$term=="gender_ratio"]<-19
model6$order[model6$term=="gender_ratio"]<-19

model2$order[model2$term=="pct_students"]<-20
model4$order[model4$term=="pct_students"]<-20
model6$order[model6$term=="pct_students"]<-20

#model2$order[model2$term=="uni_pop"]<-21
#model4$order[model4$term=="uni_pop"]<-21
#model6$order[model6$term=="uni_pop"]<-21


model1$estimate<-round(as.numeric(model1$estimate), 3)
model2$estimate<-round(as.numeric(model2$estimate), 3)
model3$estimate<-round(as.numeric(model3$estimate), 3)
model4$estimate<-round(as.numeric(model4$estimate), 3)
model5$estimate<-round(as.numeric(model5$estimate), 3)
model6$estimate<-round(as.numeric(model6$estimate), 3)


model1$std_er_star<-paste("(", round(as.numeric(model1$rse),3), ")",
                          ifelse(as.numeric(model1$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model1$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model1$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model2$std_er_star<-paste("(", round(as.numeric(model2$rse),3), ")",
                          ifelse(as.numeric(model2$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model2$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model2$pv_rse)<0.05,"*",
                                               ""))), sep = "")
model3$std_er_star<-paste("(", round(as.numeric(model3$rse),3), ")",
                          ifelse(as.numeric(model3$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model3$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model3$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model4$std_er_star<-paste("(", round(as.numeric(model4$rse),3), ")",
                          ifelse(as.numeric(model4$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model4$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model4$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model5$std_er_star<-paste("(", round(as.numeric(model5$rse),3), ")",
                          ifelse(as.numeric(model5$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model5$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model5$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model6$std_er_star<-paste("(", round(as.numeric(model6$rse),3), ")",
                          ifelse(as.numeric(model6$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model6$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model6$pv_rse)<0.05,"*",
                                               ""))), sep = "")


model1$con_std_er_star<-paste("[", round(as.numeric(model1$con_std_error),3), "]",
                              ifelse(as.numeric(model1$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model1$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model1$con_p_value)<0.05,"*",
                                                   ""))), sep = "")


model2$con_std_er_star<-paste("[", round(as.numeric(model2$con_std_error),3), "]",
                              ifelse(as.numeric(model2$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model2$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model2$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model3$con_std_er_star<-paste("[", round(as.numeric(model3$con_std_error),3), "]",
                              ifelse(as.numeric(model3$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model3$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model3$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model4$con_std_er_star<-paste("[", round(as.numeric(model4$con_std_error),3), "]",
                              ifelse(as.numeric(model4$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model4$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model4$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model5$con_std_er_star<-paste("[", round(as.numeric(model5$con_std_error),3), "]",
                              ifelse(as.numeric(model5$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model5$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model5$con_p_value)<0.05,"*",
                                                   ""))), sep = "")

model6$con_std_er_star<-paste("[", round(as.numeric(model6$con_std_error),3), "]",
                              ifelse(as.numeric(model6$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model6$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model6$con_p_value)<0.05,"*",
                                                   ""))), sep = "")


mode1_b<-subset(model1, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_ver)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_ver)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_ver)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_ver)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "Yes", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_ver)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo1_ver), "", "", "")

mode1_b<-rbind(mode1_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode1_b$model<-"Model1"
mode1_b$order<-as.numeric(mode1_b$order)
names(mode1_b)<-paste0(names(mode1_b), "_1")


mode2_b<-subset(model2, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_ver)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_ver)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_ver)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_ver)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_ver)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo8_ver), "", "", "")

mode2_b<-rbind(mode2_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode2_b$model<-"Model2"
mode2_b$order<-as.numeric(mode2_b$order)
names(mode2_b)<-paste0(names(mode2_b), "_2")


mode3_b<-subset(model3, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_brid)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_brid)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_brid)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_brid)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "Yes", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_brid)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo1_brid), "", "", "")

mode3_b<-rbind(mode3_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode3_b$model<-"Model3"
mode3_b$order<-as.numeric(mode3_b$order)
names(mode3_b)<-paste0(names(mode3_b), "_3")


mode4_b<-subset(model4, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_brid)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_brid)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_brid)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_brid)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_brid)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo8_brid), "", "", "")

mode4_b<-rbind(mode4_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode4_b$model<-"Model4"
mode4_b$order<-as.numeric(mode4_b$order)
names(mode4_b)<-paste0(names(mode4_b), "_4")


mode5_b<-subset(model5, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_bond)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_bond)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_bond)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_bond)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "Yes", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_bond)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo1_bond), "", "", "")

mode5_b<-rbind(mode5_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode5_b$model<-"Model5"
mode5_b$order<-as.numeric(mode5_b$order)
names(mode5_b)<-paste0(names(mode5_b), "_5")

mode6_b<-subset(model6, select = c("order", "term", "estimate", "std_er_star", "con_std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_bond)[1], 3), "", "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_bond)[2], 3), "", "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_bond)[3], 3), "", "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_bond)[1], 3), "", "", "")
county_fe<-c(27, "county_fe", "No", "", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_bond)$adj.r.squared, 3), "", "", "")
obs<-c(31, "obs", nobs(mo8_bond), "", "", "")

mode6_b<-rbind(mode6_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode6_b$model<-"Model6"
mode6_b$order<-as.numeric(mode6_b$order)
names(mode6_b)<-paste0(names(mode6_b), "_6")


#Step1: Collecting the models
list_order<-c(mode1_b$order_1, mode2_b$order_2, mode3_b$order_3,
              mode4_b$order_4, mode5_b$order_5, mode6_b$order_6)
#Step2: Ordering rows
list_order<-unique(sort(list_order))
list_order_df<-data.frame(list_order)
names(list_order_df)<-"order"
list_order_df$order<-as.numeric(list_order_df$order)


new_df2<-left_join(list_order_df, as.data.frame(mode1_b), by = c("order"="order_1"))


new_df3<-left_join(new_df2, as.data.frame(mode2_b), by = c("order"="order_2"))
new_df3$model<-new_df3$model_1
new_df3$model<-ifelse(is.na(new_df3$model), new_df3$model_2, new_df3$model)
new_df3$term<-new_df3$term_1
new_df3$term<-ifelse(is.na(new_df3$term_1), new_df3$term_2, new_df3$term)
new_df3<-subset(new_df3, select = -c(term_1, model_1, term_2, model_2))

new_df4<-left_join(new_df3, as.data.frame(mode3_b), by = c("order"="order_3"))
new_df4$model<-ifelse(is.na(new_df4$model), new_df4$model_3, new_df4$model)
new_df4$term<-ifelse(is.na(new_df4$term), new_df4$term_3, new_df4$term)
new_df4<-subset(new_df4, select = -c(term_3, model_3))

new_df5<-left_join(new_df4, as.data.frame(mode4_b), by = c("order"="order_4"))
new_df5$model<-ifelse(is.na(new_df5$model), new_df5$model_4, new_df5$model)
new_df5$term<-ifelse(is.na(new_df5$term), new_df5$term_4, new_df5$term)
new_df5<-subset(new_df5, select = -c(term_4, model_4))

new_df6<-left_join(new_df5, as.data.frame(mode5_b), by = c("order"="order_5"))
new_df6$model<-ifelse(is.na(new_df6$model), new_df6$model_5, new_df6$model)
new_df6$term<-ifelse(is.na(new_df6$term), new_df6$term_5, new_df6$term)
new_df6<-subset(new_df6, select = -c(term_5, model_5))

new_df7<-left_join(new_df6, as.data.frame(mode6_b), by = c("order"="order_6"))
new_df7$model<-ifelse(is.na(new_df7$model), new_df7$model_6, new_df7$model)
new_df7$term<-ifelse(is.na(new_df7$term), new_df7$term_6, new_df7$term)
new_df7<-subset(new_df7, select = -c(term_6, model_6))


new_df7$formal<-NA
new_df7$formal[new_df7$term=="vereine_tot_pop_dummy"]<-"All Networks"
new_df7$formal[new_df7$term=="bridging_tot_pop_dummy"]<-"Bridging"
new_df7$formal[new_df7$term=="bonding_tot_pop_dummy"]<-"Bonding"
new_df7$formal[new_df7$term=="post_treat"]<-"Post Week 11 Dummy"
new_df7$formal[new_df7$term=="post_treat_vereine_tot_pop_dummy"]<-"Post Week 11 x All Networks"
new_df7$formal[new_df7$term=="post_treat_bridging_tot_pop_dummy"]<-"Post Week 11 x Bridging"
new_df7$formal[new_df7$term=="post_treat_bonding_tot_pop_dummy"]<-"Post Week 11 x Bonding"
new_df7$formal[new_df7$term=="lag_covid_pc"]<-"Lag Covid per Cap."
new_df7$formal[new_df7$term=="pct_turnout"]<-"Pct. Turnout"
new_df7$formal[new_df7$term=="log_pop_total"]<-"Log Pop. Total"
new_df7$formal[new_df7$term=="log_gdp_per_cap"]<-"Log GDP/capita"
new_df7$formal[new_df7$term=="pop_den"]<-"Pop. Density"
new_df7$formal[new_df7$term=="pct_college"]<-"Pct. College"
new_df7$formal[new_df7$term=="pct_pop_above_60"]<-"Pct. Pop. above 60"
new_df7$formal[new_df7$term=="pct_pop_under_35"]<-"Pct. Pop. below 35"
new_df7$formal[new_df7$term=="pct_service_sec_narrow"]<-"Pct. Empl. Hospitality and Transport"
new_df7$formal[new_df7$term=="pct_manufacturing"]<-"Pct. Empl. Manufacturing"
new_df7$formal[new_df7$term=="east_ger"]<-"East Germany"
new_df7$formal[new_df7$term=="gender_ratio"]<-"Gender Ratio"
new_df7$formal[new_df7$term=="pct_students"]<-"Pct. Students"
#new_df7$formal[new_df7$term=="uni_pop"]<-"University-Population Ratio"
new_df7$formal[new_df7$term=="mean_mob"]<-"Mean Mobility"
new_df7$formal[new_df7$term=="sd_mob"]<-"SD Mobility"
new_df7$formal[new_df7$term=="cross_sections"]<-"No. Cross-Sections"
new_df7$formal[new_df7$term=="time_units"]<-"No. Time Units"
new_df7$formal[new_df7$term=="county_fe"]<-"County FE"
new_df7$formal[new_df7$term=="week_fe"]<-"Week FE"
new_df7$formal[new_df7$term=="weekxland_fe"]<-"Week X Land FE"
new_df7$formal[new_df7$term=="adj_rsq"]<-"Adj. R sq."
new_df7$formal[new_df7$term=="obs"]<-"Observations"



new_df8<-subset(new_df7, select = c(formal, estimate_1, std_er_star_1, con_std_er_star_1,
                                    estimate_2, std_er_star_2, con_std_er_star_2,
                                    estimate_3, std_er_star_3, con_std_er_star_3,
                                    estimate_4, std_er_star_4, con_std_er_star_4,
                                    estimate_5, std_er_star_5, con_std_er_star_5,
                                    estimate_6, std_er_star_6, con_std_er_star_6))


#Step3 Order columns
final1<-new_df8
new_row<-c("Variable", rep(c("Estimate", "Robust SE", "Conley SE"), 6))
final2<-rbind(new_row, final1)


#Step7 centering coefficients
object_latex<-xtable(final2, type = "latex", 
                     caption = "Results",
                     #Note there are X columns to be aligned below.
                     #This is because of the index column.
                     align = "llllllllllllllllllll")


comment <- list(pos = list(0), command = NULL)
comment$pos[[1]] <- c(nrow(object_latex))
comment$command <- c(paste("\\bottomrule\n",
                           "\\multicolumn{19}{l} {\\parbox[t]{47cm}\n",
                           "{\\footnotesize \\textit{Note}: 
This table shows the effect of networks on mobility in the post-lockdown period. 
The dependent variable is mobility as measured by the government. 
The interactions between the period after week 11 and different types of networks (all networks, bridging, bonding)
are the variables of interest. All networks, bridging or bonding take the value of 1 for counties with 
density of these individual networks above the mean. Similarly, the post-lockdown period is marked with 1. 
The standard errors are clustered at a county level in round brackets and 
adjusted for spatial (all networks - 455 km; bridging - 538 km; bonding - 455 km)
and serial (10 weeks) correlation in the square brackets. Significant at ∗10\\%, ∗∗5\\%, and ∗∗∗1\\%.}} \\\\ \n", sep = ""))



#Step8 Printing latex tex
tableLines<-print(object_latex,
                  #The following line controls where the lines are drawn
                  hline.after=c(-1, 1, 21),
                  floating = FALSE,
                  file = "./paper/tables/table_A4.tex", 
                  #caption.placement = "top",
                  add.to.row = comment,
                  booktabs = TRUE,
                  include.colnames=F,
                  include.rownames=FALSE,
                  sanitize.text.function=function(x){x})

multicolumns <- "& \\\\multicolumn{3}{c}{Model 1} & \\\\multicolumn{3}{c}{Model 2} & 
\\\\multicolumn{3}{c}{Model 3} & \\\\multicolumn{3}{c}{Model 4} & \\\\multicolumn{3}{c}{Model 5}
& \\\\multicolumn{3}{c}{Model 6}\\\\\\\\"
clines<-"\\\\cmidrule(lr){2-4} \\\\cmidrule(lr){5-7} \\\\cmidrule(lr){8-10} \\\\cmidrule(lr){11-13}
\\\\cmidrule(lr){14-16} \\\\cmidrule(lr){17-19}\\\\\\\\"
tableLines <- sub ("\\\\toprule\\n", paste0 ("\\\\toprule\n", multicolumns, clines, "\n"), tableLines) ## booktabs = TRUE
#tableLines <- sub ("\\\\hline\\n",   paste0 ("\\\\hline\n",   multicolumns, "\n"), tableLines) ## booktabs = FALSE
writeLines(tableLines, con = "./paper/tables/table_A4.tex", )

#SIMPLE TABLE
to_keep = c("Post Week 11 x All Networks", "Post Week 11 x Bridging", "Post Week 11 x Bonding",
            "Mean Mobility", "SD Mobility", "No. Cross-Sections", "No. Time Units", "County FE",
            "Week FE", "Week X Land FE", "Adj. R sq.", "Observations")
new_df8_short<-subset(new_df8, formal %in% to_keep)
new_df8_short <- insertRows(new_df8_short, 2 , new = NA)
new_df8_short <- insertRows(new_df8_short, 3 , new = NA)
new_df8_short <- insertRows(new_df8_short, 5 , new = NA)
new_df8_short <- insertRows(new_df8_short, 6 , new = NA)
new_df8_short <- insertRows(new_df8_short, 8 , new = NA)
new_df8_short <- insertRows(new_df8_short, 9 , new = NA)
#Taking RSE and putting in the right row
new_df8_short[2,2]<-new_df8_short[1,3]
new_df8_short[3,2]<-new_df8_short[1,4]
new_df8_short[2,5]<-new_df8_short[1,6]
new_df8_short[3,5]<-new_df8_short[1,7]

new_df8_short[5,8]<-new_df8_short[4,9]
new_df8_short[6,8]<-new_df8_short[4,10]
new_df8_short[5,11]<-new_df8_short[4,12]
new_df8_short[6,11]<-new_df8_short[4,13]

new_df8_short[8,14]<-new_df8_short[7,15]
new_df8_short[9,14]<-new_df8_short[7,16]
new_df8_short[8,17]<-new_df8_short[7,18]
new_df8_short[9,17]<-new_df8_short[7,19]

col_to_keep = c("formal", "estimate_1", "estimate_2", 
                "estimate_3", "estimate_4", "estimate_5", "estimate_6")
new_df8_short2<-subset(new_df8_short, select =col_to_keep)

new_row<-c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")
new_df8_short2 <- insertRows(new_df8_short2, 17 , new = new_row)
new_row2<-c("", "(1)", "(2)", "(3)", "(4)", "(5)", "(6)")
new_df8_short3 <- insertRows(new_df8_short2, 1 , new = new_row2)


#Step7 centering coefficients
object_latex<-xtable(new_df8_short3, type = "latex", 
                     caption = "Results",
                     #Note there are X columns to be aligned below.
                     #This is because of the index column.
                     align = "llllllll")

comment <- list(pos = list(0), command = NULL)
comment$pos[[1]] <- c(nrow(object_latex))


comment$command <- c(paste("\\bottomrule\n",
                           "\\multicolumn{7}{l} {\\parbox[t]{18cm}\n",
                           "{\\footnotesize \\textit{Note}: 
This table shows the effect of networks on mobility in the post-lockdown period. 
The standard errors are clustered at a county level in round brackets and 
adjusted for spatial (all networks - 455 km; bridging - 538 km; bonding - 455 km)
and serial (10 weeks) correlation in the square brackets.
Significant at *10\\%, **\\%, and ***1\\%.}} \\\\ \n", sep = ""))




#Step8 Printing latex tex
print(object_latex,
      #The following line controls where the lines are drawn
      hline.after=c(-1, 1, 10),
      floating = FALSE,
      file = "./paper/tables/table_1.tex", 
      #caption.placement = "top",
      add.to.row = comment,
      booktabs = TRUE,
      include.colnames=F,
      include.rownames=FALSE,
      sanitize.text.function=function(x){x})


